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Turbulent mixing and entrainment at the boundary of a cloud is studied by means of direct nu- 
merical simulations that couple the Eulerian description of the turbulent velocity and water vapor 
fields with a Lagrangian ensemble of cloud water droplets that can grow and shrink by condensa- 
tion and evaporation, respectively. The focus is on detailed analysis of the relaxation process of 
£Nj ■ the droplet ensemble during the entrainment of subsaturated air, in particular the dependence on 

turbulence time scales, droplet number density, initial droplet radius and particle inertia. We find 
that the droplet evolution during the entrainment process is captured best by a phase relaxation 
04 ■ time that is based on the droplet number density with respect to the entire simulation domain 

and the initial droplet radius. Even under conditions favoring homogeneous mixing, the probability 
density function of supersaturation at droplet locations exhibits initially strong negative skewness, 
consistent with droplets near the cloud boundary being suddenly mixed into clear air, but rapidly 
£N| 1 approaches a narrower, symmetric shape. The droplet size distribution, which is initialized as per- 

fectly monodisperse, broadens and also becomes somewhat negatively skewed. Particle inertia and 
gravitational settling lead to a more rapid initial evaporation, but ultimately only to slight deple- 
tion of both tails of the droplet size distribution. The Reynolds number dependence of the mixing 
process remained weak over the parameter range studied, most probably due to the fact that the 
inhomogeneous mixing regime could not be fully accessed when phase relaxation times based on 
H ■ global number density are considered. 

O 

^ I. INTRODUCTION 

Mixing of a passive scalar in a turbulent flow is an archetypical problem in the study of turbulence and intermittency. 
The mixing of cloudy and clear air adds additional complexity, not only because the scalar fields (e.g., droplet number 
. density, water vapor concentration, and temperature) are no longer 'passive' because of latent heating effects, but 
J> ■ also because the condensed and vapor phases are coupled through mass conservation, and because the condensed 
£NJ , phase itself can respond through various pathways. For example, upon mixing of cloudy and clear air, droplets will 
■ evaporate until the mixture becomes saturated (assuming the initial presence of sufficient condensed water), but this 
could occur by all droplets evaporating by the same amount, or by a subset of droplets evaporating completely, leaving 
the remaining droplets unchanged. The consequences for cloud properties are significant: droplet collision rates and 
cloud optical properties depend strongly on the shape of the droplet size distribution and the droplet number density. 

The entrainment of clear air and its mixing with cloudy air occurs during the entire life of a cloud. It introduces 
strong inhomogeneities at spatial scales ranging from 10 2 m down to 1mm and at time scales from hours to seconds 1 . 
Through dilution, evaporation, and perhaps enhanced collision, the entrainment and mixing process changes the 
water droplet size distribution, which is of direct consequence to rain formation and cloud radiative properties and will 
eventually determine the cloud lifetime itself. At stratocumulus cloud top, for example, entrainment and evaporation 
influence the entire cloud dynamics 3 . Here, the presence of strong wind shear can additionally enhance the entrainment 
rate^. The mixing of clear and cloudy air can be characterized by the Damkohler number, the ratio of a fluid time 
scale to a characteristic thermodynamic time scale associated with the evaporation process (phase relaxation time) 

Da=^^. (1) 

Tphase 

The two limits, Da <C 1 and Da 3> 1 characterize the homogeneous and inhomogeneous mixing, respectively. This 
notion was inspired by early laboratory experiments and the analogy with reacting flows&~— : Homogeneous mixing 
occurs when the condensational growth or evaporation of cloud water droplets is slow compared to the mixing time, 
and therefore takes place in a well-mixed environment; Inhomogeneous mixing occurs when the evaporation proceeds 
much faster than the flow structures evolve. Both processes can coexist in a turbulent cloud since a whole spectrum 
of fluid time scales is present^. 

The aim of the present work is to gain a deeper understanding of the initial evolution of the mixing processes in a 
small subvolume at the clear air-cloud interface and to characterize the multiple-scale mixing processes. Specifically, 
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we investigate how the droplet size distribution evolves as the mixing progresses. We conduct therefore a series of 
three-dimensional (3D) direct numerical simulations (DNS) which combine the Eulerian description of continuum 
fields such as velocity and vapor content with the Lagrangian evolution of an ensemble of cloud water droplets. In 
particular, we want to study the relaxation of the condensed phase during the entrainment process as a function of 
the initial droplet radius and the number density of the droplets. These two variables determine the phase relaxation 
time scale, and therefore influence the relative length and time scales at which homogeneous and inhomogeneous 
mixing predominate. Furthermore, we study the impact of turbulence on the entrainment and mixing by conducting 
simulations at different Reynolds number for the same initial vapor and droplet configuration. Finally, the effect of 
droplet inertia and gravitational settling is investigated. In order to focus on the one-way response of the droplet field 
to the turbulent mixing, we study the mixing processes in a simplified setting of model equations: the temperature is 
held fixed at a reference value of T in the simulation domain, thus leaving the saturation vapor mixing ratio constant. 
As a consequence, the turbulence is not driven by buoyancy effects as in similar studies^—, but by a volume forcing 
that mimics a cascade of kinetic energy from larger scales obeying realistic amplitudes of the turbulent fluctuations 
that match the field measurements with ACTOS platform in 9 . The reason for this choice is to disentangle the role 
of different processes on the cloud water droplet dynamics. The focus on Lagrangian behavior extends the prior 
computational studies of cloud mixing^ ' n i 12 , allowing questions of variability in droplet growth history and droplet 
inertia to be directly addressed. The Lagrangian perspective has been taken in other cloud studie d 10 ' 13 , primarily with 
an emphasis on the bulk dynamics in a turbulent cloud and the resulting supersaturation field and droplet response 
(e.g., advection-diffusion equation for the supersaturation field). This study is focused rather on the microphysical 
response to a transient mixing event. We continue in the following subsections by considering the time scales that are 
thought to govern the nature of that transient response. 

a. Fluid time scale Turbulent flows are characterized by a continuous range of time scales that can be associated 
with differently sized vortex structures or shear layers present in the flow. The largest time scale is the large scale 
eddy turnover time T — Li nt /u rms where Li nt denotes a characteristic large (energy injection) scale of the flow and 
u rms is the root-mean square of the turbulent velocity fluctuations. The smallest mean time scale is the Kolmogorov 
time Tjj = y/vjje) with the kinematic viscosity v and the mean kinetic energy dissipation rate (e) . For constant T p hase 
the spectrum of possible Damkohlcr numbers thus spans a range 

Da v = — ^— « Da < Da L = . (2) 

Tphase tphase 

The crossover from homogeneous to inhomogeneous mixing is expected to be present at Da ~ I and this can be 
associated with a length scale in the turbulent cloud^. Together with ti = Ijvg = ^ 2 / 3 /(e) 1 / 3 one gets 



Da~l t c ~^{e)T* hase . (3) 

In our simulations with constant (e) it therefore follows that the phase relaxation determines the spatial transition 
from inhomogeneous mixing (Da 3> 1) at large scales, to homogeneous mixing (Da <C 1) at small scales. Looked at 
from a somewhat different perspective, we find that the transition scale l c relative to the Kolmogorov length scale is 
simply related to a power of Da n : 

^ ~ ( V a ~ Da' 3 / 2 . (4) 

b. Phase relaxation time scale The phase relaxation time is the exponential time scale associated with the con- 
densational growth or evaporation of a population of droplets 1 ^. We provide a derivation here because it is not 
necessarily familiar within the turbulence community, yet it is of central importance in the cloud mixing problem. We 
have tried to simplify the derivation sufficiently that the key assumptions are explicitly stated, but are not obscured by 
unnecessary details; in this regard the reader is also referred to Kostinski's lucid treatment^. We begin by considering 
a single cloud droplet, assumed to be in equilibrium with its surrounding vapor field, so a steady mass flux of vapor 
toward the droplet surface is accompanied by a steady flux of latent heat away from the droplet surfac o 14 ' 16 . Both 
steady fluxes are given by Fick's law 

Fv = _ D ^IL F = _ fc — (5) 
dR ' Q dR ' w 

with the vapor mass density p v , the mass diffusivity D, the temperature T, and the thermal conductivity k. At the 
droplet radius r, we define boundary conditions T(R = r) = T r and p v (R = r) — p v ^ r . In steady state the resulting 
profiles reach their asymptotic values of a reference temperature T(R 3> r) — and a reference vapor mass density 
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p v (R ^> r) = p V: oo at several droplet radii. Later it will be assumed that the reference values are the same for all 
droplets in a simulation 'grid box,' which is an implicit statement that fluctuations in the vapor and temperature 
fields due to 'local' droplet interactions are neglected^. A steady regime requires energy conservation, i.e., 

F v + LF Q = , (6) 

with L being the latent heat of vaporization (see Table 1). Inserting Eq. ([5]) into balance (|6]) and integrating with 
respect to R results in a relation between mass density values and temperatures 

Pv,co pv,r k 

T r Too LD 

It is further assumed that the water vapor pressure at the droplet surface is at the saturation value e s (T r ) and thus 
Pv,r — Pvs ■ Both are connected via the ideal gas law e s — R v p vs T where R v is the vapor gas constant given in Table 
1. By translating the linearized solution of the saturation pressure from the Clausius-Clapeyron equation into an 
expression for the saturation vapor mass density, ones arrives together with T r s» and (J7]) at 

frp \ (T \ r*j LpvsjTco) {rr _rp \ _ LPvJ^Oo) LP . . , . 

Pvs\-loo) Pvs\J-r) — rp 2 \ °° r > ~ f> rp2 \P v , r Pv,oo) ■ (°) 

The vapor flux across a sphere of radius R has to be equal to the change of liquid water mass inside the sphere 

„q „ do,, dMi ,dr 



Integration from R = r to R = oo results to 



dr D 

^37 = — \Pv,oo - Pv,r) ■ (10) 
at pi 



With Eq. (|8]), we can substitute p v ,r(= Pvs{T r )) in Eq. dTUl) and get the following equation for the radius growth by 
condensation 

dr D ( DtfpnjT^ y 1 V 

r ~n- — {Pv,oo - Pvs{Too)) [ 1 -\ , p ^ ^ — [Pv,oo ~ Pv S [Too)) ■ (11) 

at pi \ kRvlSo J Pi 

The diffusivity constant T> now incorporates the self-limiting effects of latent heat release This modified diffusivity 
can be written in terms of a ratio of two heat fluxes T>/D = (1 + ^/^fc) -1 : A characteristic heat flux due to 
latent heating resulting from a small change in droplet temperature, $l = LDAp vs , where A.p vs is obtained from 
AT through the linearized Clausius-Clapeyron equation; And a heat flux due to thermal conduction for the same 
temperature difference, $£. = kAT. For typical warm cloud conditions $L/$k is of order unity, so the heat-transfer- 
limited diffusivity T> can be reduced by a factor of 2 (as is the case for the example values in Table 1). Only at rather 
low temperatures, e.g. < —20 °C, do the thermal effects become negligible for liquid water. 

The phase relaxation time scale is a direct consequence of the combination of Eqns. ^ and (|TT|) 

— — = 4rfDr(p v>ao -pvsiTca)), (12) 
at 

and the conservation of water mass, n^dAf/dt = — dp VtOC /dt. Here rid is the droplet number density and thus 

dpv,oo 



dt 



-ATmaVripy.^ - p vs (T OQ )) . (13) 



If Too is constant, such that p 1)s (T 00 ) is constant, this equation describes an exponential relaxation with a characteristic 
time constant, the phase relaxation time, given by 

Tphase = 4nn d T>r ' (14) 

It should be noted that a more detailed derivation expresses the phase relaxation time in terms of the first moment 
of the droplet size distribution (integral radius). Furthermore, we specifically call attention to the assumption of 
uniform number density, which of course is not exactly the case during an inhomogeneous mixing event in which 
the local number density varies considerably. This is one of the motivations for taking a Lagrangian perspective, 
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where individual droplets are followed through the flow, as opposed to treating the condensed phase as a continuous 
mediu m 11 1 12 . 

The outline of the manuscript is as follows. After introducing the set of equations which are solved numerically in the 
present Euler-Lagrangian model, we will discuss in brief some properties of the statistically stationary turbulent state. 
This discussion is followed by a description of the initial vapor content profile. The cntrainmcnt process is described 
afterwards in combination with an analysis of the mean volume radius, the size distribution of the droplets and the 
supersaturation along the droplet trajectories. Finally, we will add the effect of particle inertia and gravitational 
settling to the simulations and quantify their impact. We conclude the work with a summary and an outlook. 



II. MODEL EQUATIONS AND NUMERICAL METHOD 



The turbulent velocity field u(x, t) and the pressure field p(x, t) are those necessary for the description of an 
incompressible turbulent flow. In this flow the vapor mixing ratio field <z„(x, t) is transported and diffuses. The vapor 
mixing ratio is defined as 

q v {x,t) = — , (15) 
Pd 

where p v and pd are the mass densities of vapor and dry air, respectively. For the purposes of this paper, the 
advection-diffusion equation for temperature is not considered. 
The Eulerian equations for the turbulent fields are 

Vu = 0, (16) 

<9 t u+(u-V)u = — Vp + iA7 2 u + f, (17) 
Po 

d t q v +u-\7q v = DV 2 q v ~C dl (18) 

where f (x, t) is a bulk forcing which sustains the turbulence and Cd is the condensation rate. The entrainment is stud- 
ied in a cube with volume V = l? x and with periodic boundary conditions in all three spatial directions. It is spanned 
by an equidistant mesh with iVj cells of mesh size a. The Eulerian equations are solved by a pseudospectral method 
using fast Fourier transformations. Time advancement is done by a second-order predictor-corrector method. The 
spectral resolution in the present cases is k max r\ = 3 with the maximum resolved wavenumber k max — 2tt^/2N x /(3L x ) 
and the Kolmogorov scale r\ (see Table 1). Grid sizes used throughout this work are N% = 128 3 , 256 3 and 512 3 
corresponding with turbulent flows at Taylor microscale Reynolds numbers R\ = 42, 59 and 89, respectively. 

c. Volume forcing Here, we consider a turbulent flow that is sustained by a volume forcing f (x, t) in a statistically 
stationary turbulent state. This driving is implemented in the Fourier space for some modes with the smallest 
wavenumbers fc/ only, i.e. kj 1 ss L x . The kinetic energy is injected at a fixed rate €i n into the flow. The volume 
forcing is established by the expression^ 

f(M) = ei W = U( , k '.? T^ ftc.k,, (19) 



with the Kronccker delta 



Jk.k f 



•£k, e x>(k/>*)l 5 



1 if k = k/, 5k. k/ = otherwise, (20) 



and the wavevector subset K. which contains some wave vectors, e.g. k/ = (1, 1, 2) plus all permutations with respect 
to components and signs. Since the large-scale velocity follows a Gaussian statistics, the present forcing will act in 
similar way as other stochastic forcing schemes^. The present energy injection mechanism prescribes the mean energy 
dissipation rate; that is, the magnitude of the first moment of the energy dissipation rate field, (e), is determined 
by the injection rate, e;„, having no Reynolds number dependence. This can be seen as follows. Given the periodic 
boundary conditions in our system, the turbulent kinetic energy balance, which results from rewriting (|17p in the 
Fourier space, follows to 

= -,£fc 2 |u(k,i)| 2 +£f(k,t).u*(M), (21) 
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TABLE I: List of constants and reference values are given in the upper part of the table. Simulation parameters and charac- 
teristics of the statistically stationary turbulent state follow in the second one. 



Quantity 


Symbol 


Unit 


Value 


Reference temperature 


T 

' oo 


K 


270 


Reference pressure 


Poo 


hPa 


845 


Kinematic viscosity 


V 


2 — 1 

ma 


1.5 x 10~ 5 


Vapor diffusivity at Too 


D 


9 — 1 

m 2 s 1 


2.16 x 10~ 5 


Modified vapor diffusivity at Too 


V 


2 —1 

m s 


1.31 x 10~ 5 


Thermal conductivity of air at Too 


k 


t — 1 -lr/- 

Jm a K 


1 2.38 x 10~ 2 


Gravity acceleration 


a 


m s 


9.81 


Gas constant for water vapor 


R v 


J K kg 


461.5 


Gas constant for dry air 


Rd 


JK kg 


287.0 


Specific heat at constant pressure 


Op 


Jkg^K' 1 


1005 


Latent heat 


L 


J kg- 1 


2.5 x 10 6 


Liquid water density 


Pi 


7 — 3 

kgm 


l(r 


Reference mass density of air 


PO 


kgm 


1.06 


Saturation pressure at Too 


e s (Too) 


Pa 


484 


Saturation vapor density at Too 


pvs (Too ) 


kg m~ 3 


q n w i n — 3 

o.y x tu 


Constant in Eq. ([26]) 


K 


2 —1 

m s 


5.07 x 10 


Box length 




m 


U.lzo, U.zOD, U.olz 


Grid resolution 




mm 


1.0 


Kolmogorov scale 




mm 


1.0 


Mean energy dissipation rate 


(e> (= e ln ) 


m 2 s~ 3 


0.003375 


Root-mean-square velocity 




cm s _1 


8.6, 10.1, 12.5 


Taylor microscale Reynolds number R\ = ^5/(3v(e)) u%. ms 




42, 59, 89 


Initial cloud water droplet radius 


Ro = r(t = 0) 


fim 


10, 15, 20 


Cloud number density 


rid 


cm' 3 


62, 82, 164, 328 



where u* is the conjugate complex Fourier mode. The first term on the right hand side of (|21j) is the volume average 
of the energy dissipation rate. Additional time averaging in combination with (I19[) results in 



vJ2k 2 (Hk,t)\ 2 ) t = (e) =e in = £(f(k,t).u*(M)) t . 

k 



(22) 



The applied driving thus allows a full control of the mean energy dissipation rate (e) and thus the Kolmogorov scale 
r\ via the parameter e in in (fT9|) . 

d. Cloud water droplet advection and condensation rate The Lagrangian evolution of each of the N droplets in 
the volume V is described by the following set of equations 



dX(i) 
dt 

dt 



= v(<), 

1 



[u(X,t)-V(t)]+g. 



(23) 
(24) 



Here, X is the droplet position and V its velocity. We consider both droplets which match perfectly with the 
surrounding fluid velocity as well as inertial particles with a finite particle response time t p = Ipir 1 /(9po,v). 

As droplets are advected by the fluid they can grow or evaporate in response to the local vapor field (recalling that 
in this study temperature is constant). Direct droplet interactions through collision are neglected in order to focus 
solely on the initial stage of the entrainment and mixing process. The vapor mixing ratio can be coupled to droplet 
growth by defining the supersaturation S — p v ,oc/ Pvs(Too) — 1, such that 



S(x,t) 



Qv,. 



(25) 



Then it follows from Eqn. (fTTj) in Sec. Uthat the droplet growth rate can be written as 



dr 

r— = KS 
dt 



with 



K 



Pi 



L 2 



\De s {T 00 ) kRyT^ 



(Too) 

Pi 



V. 



(26) 
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FIG. 1: Turbulent kinetic energy (left) and volume averaged kinetic energy dissipation rate (right) as a function of time. The 
dashed line in the right panel marks the prescribed ensemble average {s)v,t (see Table 1). For t > 15s, both quantities and 
the turbulence as a whole are fully relaxed into a statistically stationary state. This is the starting point of the entrainment 
simulation which is marked by a vertical dotted line. Data are for the run at R\ = 89. 



In the Lagrangian frame this condensational growth process becomes 



r(t)^ = KS(K,t). (27) 



We calculate the condensation rate field C<j(x, t) following— by 

A 

K 

dt poa 



r< ■ m 1 drrLl ( x > f ) kxpiK v 

c*j(M) = - ^ — = ^Tl^SOtp^Xt)' (28) 

p=l 

where m a is the mass of air per grid cell and the sum collects the droplets inside each of the grid cells of size a 3 
that surround the (grid) point x. This relation closes the system of Eulerian-Lagrangian equations. The transmission 
of the Eulerian field values at grid positions to the enclosed droplet position is done by trilinear interpolation. The 
inverse procedure is required for the calculation of the condensation rate which is evaluated at first at the droplet 
position and then redistributed to the nearest eight grid vertices. 

The particle advection poses a numerical problem when particle inertia is considered. Evaporating droplets cause a 
particle response time t p — > thus making Eqn. (|24[) stiff. Therefore a semi-implicit second-order particle advection 
scheme is chosen. While the equations for the radius and the droplet position are solved by a predictor-corrector 
scheme, the droplet velocity equations are solved by a combination of an implicit forward Euler step that is required 
for the corrector step of the droplet positions and a trapezoidal scheme for the velocity itself. We verified the accuracy 
of this scheme by the analytical test case of the freely falling droplet with friction for u = for Stokes numbers down 

tO Strj = Tp/Tjj r~J 10 -4 . 

The complete list of thermodynamic reference values and constants is summarized in Table 1. The values for 
reference density, temperature, pressure and the resulting saturation values are chosen in agreement with recent 
airborne measurements by Lehmann et alA It is worth emphasizing again that, as discussed in the introduction, 
we have set up the model equations such that there is no active feedback to the turbulent dynamics through the 
temperature field. Specifically, while the effects of latent heat and thermal conductivity are included in the droplet 
growth rate (e.g., through the modified diffusivity T> in Eqn. (fTTj)). there is no coupling of the condensation rate 
(which is given by Eqn. (|28l) ) to the momentum equation (1171) via a buoyancy term, and no advection-diffusion of a 
temperature field. This study is focused on the vapor advection-diffusion aspects of the mixing problem. 



III. PREPARATION OF THE INITIAL TURBULENCE STATE 



In order to prepare the turbulence initial conditions for the Euler-Lagrangian simulations we first run a pure 
flow simulation with the volume driving described by (|19p . Figure [T] demonstrates the relaxation into a statistically 
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FIG. 2: Statistics of the velocity gradients. Left: Probability density functions (PDF) of the energy dissipation rate field 
e(x, t) and the enstrophy density w 2 (x, t). Right: Comparison of the energy dissipation rate field statistics with the lognormal 
prediction of the refined similarity hypothesis 2 ^. The dotted vertical lines are at z = ±2. Data are again for the run at R\ = 89. 



stationary state by means of the time traces of the turbulent kinetic energy (left) and the volume-averaged energy 
dissipation rate (right). For times larger than 15 seconds both quantities are found to fluctuate moderately about 
their temporal means. We also checked that the isotropy of the flow is established by comparing the mean squares of 
the three velocity components. Figure [2] shows the probability density function (PDF) of the energy dissipation rate 
field and the square of vorticity magnitude, denoted as enstrophy density, which are given by 

The stretched exponential tails of both quantities demonstrate the enhanced spatial intermittency of the velocity 
gradients at the smaller scales (see e.gJ^). The right panel displays the PDF of the energy dissipation rate field (same 
data as in the left panel) in comparison with the refined similarity hypothesis prediction by Kolmogorov>2£. Deviations 
in both tails for \z\ > 2 are found. Two aspects contribute, in our view, to the deviations. First, the Reynolds number 
of the present simulations are still moderate. Second, our spectral resolution exceeds standard resolutions by at least 
a factor of 2. In Ref it was demonstrated that the higher spectral resolution is necessary to resolve the tails, i.e. the 
rare high- amplitude events, sufficiently well. We also verified from the statistical analysis that the relation (e) = v(uj 2 ) 
is satisfied. 

In Table 1 we summarize the simulation parameters and turbulence quantities. The numbers are matched to 
typical magnitudes for turbulent cumulus clouds 9 . The developed turbulent state serves as the initial condition for 
our entrainment investigations. 

By definition, mixing presumes that the initial state possesses imposed gradients in the fields, but the form of 
those gradients is by no means obvious. For example, the initial vapor-droplet configuration depends on assumptions 
regarding the degree of correlation between the two fields, as well as the spatial scales of the fluctuations. In order to 
focus on mixing from a state possessing a single, clearly defined initial length scale, we have chosen to begin with an 
idealized slab-cloud geometry. This slab could be considered analogous to a cloud edge or cloud filament boundary 
formed through inhomogeneous mixing at larger scales. More precisely, the initial condition is a slab-like filament 
of supersaturated vapor that fills about one third of the cubic simulation box, with the vapor profile across the slab 
given by 

q v (x, y,z,t = 0) = (q'r x ~ Qv) exp [~A(x ~ *o) 6 ] + ■ (30) 

Here, q™ ax is the maximum amplitude of q v , which exceeds the saturation value q vs {To) by 2%. The variable q% 
stands for the environmental vapor mixing ratio, representing the subsaturated clear air outside the supersaturated 
filament. Droplets are seeded randomly in the supersaturated slab as a monodisperse initial ensemble. Technically, a 
sharp boundary between clear and cloudy air on a length scale smaller than l c is realistic only in a transient sense, 
and other vapor profiles have been considered (see e.gJ^). Furthermore, as already stated, we use fully developed, 
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forced turbulence throughout the volume, which is idealized with respect to the turbulence level in- and outside a 
cloud (see e.gi^i). One might consider, however, that both the transient sharp boundaries and the uniformity of 
turbulence energy dissipation rate mimics the lower end of the turbulent energy cascade during a mixing event. Most 
importantly, the justification for the idealized initial conditions is that they allow us to study the entrainment in a 
well-controlled setup with clearly defined length and time scales, and therefore to disentangle the contributions of 
different physical processes to the droplet dynamics. 

We proceed now to perform a series of simulations that combine three initial radii Rq = r(t — 0) with four different 
droplet number densities rid as listed in Tab. [TXJ Corresponding to these number densities (rid), the total numbers 
(AO of advected droplets are 412.500, 550.000, 1.100.000 and 2.200.000, respectively for R x = 59. Hence, a total of 
12 simulations are run for a time interval of 60 seconds which corresponds with 60.000 simulation time steps. In 
addition, we add two runs at different Reynolds numbers. All values have been chosen so as to be representative of a 
warm cumulus cloud, but also so as to allow exploration over a realistic range of microphysical time scales (e.g., phase 
relaxation time as defined in Eqn. (|14[) ). The 'global' number density n d g ' ^ and liquid water content w( g >°\ i.e., the 
values obtained when droplet number and liquid water content are averaged over the entire computational domain, 
for all simulations are also listed in Tab. [II] Phase relaxation times based on both the initial cloud properties and the 
global averages are also listed, as well as the resulting Damkohler numbers. These will be discussed in detail later, 
but for now suffice to say that cloud values of Da n are always less than 0.1, safely in the homogeneous mixing limit, 
and cloud values of Dcll vary from approximately 0.3 to 3, thereby just making the transition toward inhomogeneous 
mixing at the large scales. Several of the liquid water contents are rather high, but this allows the largest Dcil to be 
achieved. 



IV. ENTRAINMENT PROCESS AND DROPLET EVOLUTION 



Figure [3] displays the time evolution of the turbulent entrainment process. The two isosurfaces are shown for which 
the vapor content satisfies q v (x,y, z,t) ~ q vs . Initially planar, they become highly convoluted as the time progresses. 
For times larger than 2 seconds, approximately one large eddy time, the vapor fluctuations decay below the saturation 
mixing ratio, leading to evaporation of all droplets. In Fig. 0] the same evolution is shown from the Lagrangian 
perspective. We project the position of the individual droplets of the whole ensemble onto the x-y plane. The 
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FIG. 4: (Color online) Time evolution of the droplet ensemble. From top left to bottom right: t = 0.0 s, t = 0.2 s,t = 0.7 s and 
t = 1.0 s. The positions of all 5.5 x 10 5 (which corresponds with rid = 82cm -3 ) are projected onto the x-y plane. The vertical 
lines in the first three panels indicate the minimum/maximum x position up to which vapor filaments with q = q V3 have been 
advected. Data are the same as in Fig. [3] 



additional vertical lines indicate the minimum and maximum x position up to which filaments with q v {x, y, z, t) — q vs 
have been advected by the turbulent flow. At the beginning of the evolution, droplets and vapor filaments follow the 
same stretching and twisting phenomena. With progressing time, however, both dynamics become decoupled as the 
vapor field is subjected to additional diffusion and a weak source-sink contribution from the condensation rate. This 
can be observed by comparison of the upper right (t = 0.2 s) and lower left (t = 0.7 s) panels. 

Figure [5] shows the total water balance for a typical evolution run. The sum of vapor and liquid water masses in 
the volume V ', M v and Mi, has to be constant with respect to time 

r 4 Np 

M v (t) + Mi(t) = pol q v (x,y,z,t)dV + -irpiJ^rUt) = const. (31) 
J v 6 i=l 

We can distinguish two cases for the dynamics of the present entrainment problem: 

• Case 1 is present if 



M,(i = Q) < M vs - M v (t = 0) , 



(32) 
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FIG. 5: Total water content conservation versus time, 
liquid water mass Mi. The data are for the run at R\ - 
the dashed line. 



The total mass in volume V consists of the vapor mass M v and the 
: 89. The saturation level and the corresponding M va are indicated by 



where M vs = poq vs V is the total mass of the saturated vapor. As a consequence 

Qv Qv,oo < Qvs and r t ->• Vi = 1...N . (33) 

Eventually all droplets evaporate and the vapor content fluctuations decay. The vapor content relaxes to a 
homogeneous field amplitude q v ^ below the saturation level. 

• Case 2 is present if 

Mi (t = 0) > M vs -M v (t = 0). (34) 

In this scenario the dynamics ends with the state 

Qv -> qvs and n -» r liOC Vi = l..JV. (35) 

Eventually, the vapor content fluctuations decay to zero and the vapor content relaxes to the saturation value 
q vs . Thus >S(x, t) = and all droplet growth or evaporation ceases, leaving a droplet population with a mean 
volume diameter (r 3 )^. It should further be noted that N is not necessarily conserved during the process 
because some subset of the droplet population may evaporate completely. 

In the following Lagrangian statistical analyses of the droplet growth, the supersaturation along the droplet paths 
and the relaxation time scales are discussed. Due to the initial excess of saturation level, the majority of the droplets 
can grow in the first one to two seconds. Subsequently they start to shrink slowly. The progressing entrainment of 
clear air generates local filaments of subsaturated vapor for an increasing number of droplets beginning with those 
closer to the original planar cloud-clear air interface. 

Figure [5] shows the evolution of the droplet size distribution for examples of the two scenarios just described: Rq — 
10/iTO and Rq — 15/im with droplet number density — 164cto~ 3 . The left panel of the figure illustrates the evolution 
of the size distribution for Rq = 10/zm and indicates that after t = 25s the first droplets have been evaporated. At 
this point the size distribution has developed a pronounced negative tail that appears nearly exponential. 

After a minute, the majority of the droplets have evaporated. This continues until the last droplet is vanished (not 
shown), corresponding to Case 1. The right panel of the same figure illustrates the evolution of the size distribution 
for Rq = 15/xm, corresponding to Case 2. After 60 seconds, all the droplets have stopped shrinking or growing, 
leaving a relatively narrow, but negatively-skewed size distribution. Droplets will be further advected in the sustained 
turbulent flow, but they do not shrink or grow in the homogeneously mixed vapor field. In both examples the negative 
skewness of the size distribution is primarily a result of the nonuniform exposure of droplets to subsaturated air. 

Time scales have been of central interest in the discussion of inhomogeneous versus homogeneous mixing for 
decades^— and there has been some disagreement about which time scales correctly represent the microphysi- 
cal response (e.g., see22). We address the question here directly. The steady-state saturation level in Case 2 is reached 
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TABLE II: List of parameters for the phase relaxation studies. We list the initial droplet radius Ro, the initial number densities 
with respect to the slab cloud, n^'°\ and the whole volume, n^'°\ the global liquid water content i// 9 ' - 1 in gm~ 3 , the single- 
droplet evaporation time r r , the phase relaxation times (|14[) based on the number densities n^ c ' ' and n^ 9,0 ' , the numerically 
determined relaxation time r r£ i ax , and the large-scale eddy turnover time T. Furthermore the four possible Damkohler numbers 
and the Taylor microscale Reynolds number are given. All number densities are in cm -3 and all times in s. The Kolmogorov 
time scale in all 14 DNS runs is t v = 0.067s. Empty entries for r re i ax stand for the complete evaporation of the droplets (Case 
1). 
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with an approximately exponential relaxation time r re i ax , which we compare with the phase relaxation time r p h a se as 
given by Eq. ([j~4")) and the single-droplet evaporation time r,. = —Rq/2KSo (the latter obtained via integration of Eq. 
(|26[) assuming constant S = So). For times larger than T re i ax the cubic mean droplet radius (r 3 )^, which is directly 
proportional to M/(t), becomes constant as illustrated in Fig. [7] (left). The symbol (-)l denotes an average over the 
Lagrangian droplet ensemble. The relaxation time of the droplet ensemble T re i ax is obtained by a fit of a decaying 
exponential to the simulation graphs of the mean cubic radius. The fits to the data are indicated by the symbols 
in the left panel. The various time scales for all simulations are given in Tab. [TT1 It is immediately evident that 
neither the single droplet evaporation time r r nor the cloud phase relaxation time based on the cloud values T^^' se 
correctly accounts for the observed relaxation. Figure [JJ (right) displays, instead, the comparison of the findings for 
Treiax with the phase relaxation time r^l e as calculated from (fLlj) but using the global number density n^°\ As 
shown in Fig. [7J (right), the indirect proportionality T re i ax ~ r -1 as given by the theoretical prediction is confirmed. 
The same holds for r re i ax ~ rQ . One can observe that the relaxation times T re \ ax are slightly greater than the 

corresponding values of T^'° se in most runs, presumably due to the steadily decreasing radius and the finite rate at 
which the droplets are spread throughout the volume. It should be noted, however, that calculating a T p hase based 
on the global number density and the final droplet radius was less consistent with the observed T re i ax . The simple 
phase relaxation model is therefore a surprisingly good representation of the microphysical response to the mixing 
process in the range of Damkohler numbers investigated. This range just barely approaches Da^ ~ 1, and so is most 
representative of homogeneous mixing, just approaching the transition stage in the simulation with the highest liquid 
water content. 

In order to understand the small quantitative disagreement, one has to recapitulate the assumptions that enter 
the derivation of the phase relaxation time in (|14p . There, the droplet is embedded in a homogeneous vapor field, 
an assumption that is not fully sustained in the entrainment simulation. The different Lagrangian history of each 
individual droplet and the permanently changing saturation conditions cause a slower relaxation than the idealized 
situation assumed in the derivation of the phase relaxation time. As mentioned before, the initial vapor profile 
has a maximum amplitude of the supersaturation of 2%. Due to entrainment process, the value of supersaturation 
starts decreasing from the beginning. Figure [8] (left panel) depicts the PDF of the supersaturation S(X.,t) along the 
Lagrangian cloud droplet paths, monitored at different times up to 7 seconds. The figure indicates that in the first 
seconds of the evolution the left tail of the PDF steadily grows, becoming approximately exponential. This initial 
time corresponds with the time required for clear air to reach the center of the original slab cloud, i.e., the large eddy 
time T. It is consistent with the view given by Fig. [4j in which a subset of droplets are mixed into the clear air and 
therefore experience stronger evaporation than the average. Afterwards the left tail of the PDF narrows continuously 
until all droplets have reached the same subsaturation level of about —2.5% in this particular example. Ultimately, 
the transient, nonuniform exposure during the early mixing leads to the negatively skewed size distributions shown 
in Fig. [51 still preserved long after the transients have decayed. 
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FIG. 6: Evolution of the droplet radii. The left panel represents evolution of the PDF for an initial radius Ro = 10/ira (Case 
1) and right panel for Ro = 15/im (Case 2). Both simulations started with total number of droplets of N = 1.100.000. Data 
are for = 256 3 and i? A = 59. 



V. REYNOLDS NUMBER DEPENDENCE OF THE ENTRAINMENT 



The Reynolds number dependence of the relaxation process is studied in the following way: we prepared initial 
profiles for = 62cm -3 and i?o = 20/im in volumes of side lengths L x = 12.8, 25.6 and 51.2 cm. The initial vapor 
profile pop obeys the same parameters in all three runs. Thus the volume of the initial slab cloud and the droplet 
number N increase by a factor of 8 and 64 when going from L x = 12.8 cm to L x — 25.6 cm and L x = 51.2 cm, 
respectively. 

Figure [9] indicates that the decay and thus the relaxation time differ only slightly with increasing Reynolds number. 
This holds particularly for the final phase of the relaxation. Slight differences can be observed in the initial phase 
of the entrainment. The growth of the droplets is most pronounced for the simulation in the biggest domain. The 
reason is that the ratio of entrainment area to volume is the smallest. Thus more droplets can grow unperturbed. As 
in the last section, the fit of an exponential profile to the graphs accounts for this difference and includes the part of 
the curves only in which M[(t) decays monotonically. The relaxation times obtained are r re i ax = 12.3s for R\ = 89, 
Treiax = 12.4s for R\ = 59 and r re i ax — 13.0s for R\ = 42, as given in Tab. [TTJ The relaxation becomes slightly faster 
as the Reynolds number is increased. By comparison, the global phase relaxation time for all three simulations is 
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Time in s r in urn 

FIG. 7: The phase relaxation of the droplet ensemble. Left: Mean cubic radius versus time for different initial radii and droplet 
number densities (solid lines). They are directly proportional to the changing liquid water mass in the volume. Fits to the data 
in order to determine the corresponding e-folding times are overlayed by symbols. Number densities are given in the legend. 
Right: Comparison of the theoretical prediction for the phase relaxation time T^^ e (lines) with the numerical value of r re i ax 
(symbols). Symbols agree with the left panel. All data are JVj = 256 3 and R x = 59. 

VI. ROLE OF PARTICLE INERTIA AND GRAVITATIONAL SETTLING 

The results of the droplet evolution presented in the previous section were obtained from a Lagrangian model without 
particle inertia and gravitational settling. Thus, the particle velocities are exactly the same as the velocities of the 
surrounding fluid flow, i.e. V(t) = u(X, t). In the next step, all simulations have been repeated with gravitational 
settling and a finite particle response to variations of the local advecting velocity. This results in solving the full 
set of Eqns. (|23[) (|27|) . The maximum Stokes numbers which are obtained in the simulations are of the order of 
Stv f° r = 20fim, and the corresponding settling parameter is Sv^ = St rj (g/a r) ) < 4, where a v is the 

Kolmogorov acceleration^. Therefore, we expect that effects of particle inertia and gravitational sedimentation are 
just becoming significant for the largest droplet sizes considered. 

The relatively simple dynamical change due to gravity is illustrated in the PDF of the vertical droplet velocity 
component, shown in the left panel of Fig. QTJ] Two distributions with and without inertia are depicted for i?o = 20/jto 
and rid = 164cm -3 . There is an offset to a small negative amplitude that stands for the slow downward motion of 
the droplets, but the shape of the PDF is essentially unchanged. From the right panel, we can conclude that the 
mean radius of the droplets relaxes to a value of about (r)oo ~ 19/zm which corresponds with an inertial time scale of 
T. p w 0.005s; the terminal velocity estimate is v g = —gr p sa —5 cms -1 , which is consistent with the offset observed in 
the right panel of the figure. 

But the microphysical effects of inertia and gravitational settling are more subtle. As shown in the right panel of 
Fig. [51 the initial evolution of the supersaturation PDF is significantly altered by the presence of droplet inertia and 
gravity: the left tail of the PDF grows more rapidly, presumably as a result of stronger droplet decoupling from the 
fluid, but ultimately the supersaturation PDF collapses into a similar, relatively narrow but symmetric distribution. 
Droplets are initially exposed to very different vapor environments as a result of their inertia and gravitational 
settling, but ultimately the primary influence is on the evolution of the negatively skewed tail of the supersaturation 
pdf, and the overall homogenization of the supersaturation field is unchanged. A comparison of the droplet size 
distributions further unravels the systematic difference between the droplet dynamics with and without both effects. 
Figure [10] (right) indicates a more rapid evaporation and thus a more rapid pronounced left tail of the droplet size 
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FIG. 8: Evaluation of the PDF of the supersaturation for the case of an initial radius of Rq = 20/im and a droplet number 
density of rid = 164cm -3 . Data are again for N% = 256 3 . At the beginning, supersaturation S is 2% and after 7 seconds it 
relaxes to about 2.5% subsaturation. Left: without particle inertia and gravitational settling. Right: with particle inertia and 
gravitational settling. 

distribution. At later times, both become very similar but with slight depletion in both tails. This means that 
the droplet evaporation is enhanced initially when gravitational settling and inertial effects are present, presumably 
because these effects lead to more rapid decoupling of droplets from the fluid containing high supersaturation. We 
verified, however, that the variation of the mean radius with respect to time is almost unchanged when inertia is 
included. This also can be interpreted as a result of the more rapid decorrelation of particles from the fluid, such 
that droplets that do come into contact with anomalously low supersaturations, do so for shorter times. The present 
Lagrangian approach is particularly well suited to unraveling these details of the entrainment process. To summarize 
this part, it is found that particle inertia and gravitational settling have a strong influence on the initial evolution of 
the supersaturation PDF and on the positive and negative tails of the droplet size distribution, but that the transient 
effect is sufficiently small as to have essentially no influence on the mean droplet size. 

VII. SUMMARY AND OUTLOOK 

A model for turbulent mixing and entrainment, that couples the Eulerian description of the velocity u and water 
vapor mixing ratio q v with a Lagrangian ensemble of cloud water droplets has been presented in this paper with an 
emphasis on understanding the dynamics of the turbulent entrainment at the interface between clear and cloudy air. 
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FIG. 9: Decay of the total water mass Mi(t) versus time. The curves of the two smaller Reynolds numbers have been rescaled 
by a factor of 8 and 64 in order to collapse the curves. 
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FIG. 10: Effects of particle inertia. Left: The PDF of the vertical Lagrangian velocity component at t = 10s with and without 
inertia effect. Right: Droplet size distributions taken at four different times. Solid line graphs represent size distributions with 
inertia and graphs with symbols represent the corresponding data without particle inertia as indicated in the legend. In both 
panels we took Ro = 20/im and — 164cm -3 . All data are again for TV 3 = 256 3 and times in seconds. 



The direct numerical simulation model, which resolves the turbulence in a small subvolume of the cloud down to 
the Kolmogorov length rj, is capable of generating turbulent flow conditions as observed in cumulus clouds 9 , e.g., low 
mean values of the kinetic energy dissipation rate are consistently obtained. Microphysical and turbulence parameters 
have been chosen to explore the two limiting cases of turbulent mixing in this setup, homogenous and inhomogeneous 
mixing. A central quantity to describe this physical process is the phase relaxation time, which has to be compared 
with the continuum of turbulence time scales. 

Two basic dynamical scenarios are possible in the present setup, depending on the amount of liquid water present 
at the beginning of the entrainment process. The first case leads to a complete evaporation of the droplets, and 
the second ends with a steady state droplet population surrounded by a saturated homogeneous vapor field. The 
relaxation time T re i ax to this state (which is obtained from the simulations) is compared with the phase relaxation 
time Tphase- The observed T re i ax display the expected dependencies on rid and r from Eq. (1141) . The magnitude of 
Treiax, however, is found to be significantly larger than T p h ase in all simulations, where the phase relaxation time is 
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calculated with the undiluted cloud droplet number density (as is customary in the literature, to our knowledge). In 
contrast, very close agreement with the observed r re i ax is obtained when r p /j ase is calculated using the diluted (or 
'global') number density. 

The Lagrangian approach has allowed for detailed analysis of the droplet size distribution in conjunction with the 
evolution of the supersaturation field sampled at the droplet locations (supersaturation PDF). During the transient 
mixing event the initially perfectly monodisperse droplet population broadens significantly, with a distinct negative 
skewness. This is partially a result of the strongly negatively skewed supersaturation PDF, which at early times in 
the mixing displays nearly exponential tails on the negative side of the distribution. This skewness arises from the 
droplets at the interface of the cloud, that are suddenly mixed into the clear air. Interpreted another way, this is 
an early manifestation of microphysical effects of inhomogeneous mixing, in which a subset of droplets is assumed to 
evaporate completely, leaving the remainder of droplets unchanged. These simulations have been performed primarily 
with Dcil < 1, i.e., favoring homogeneous mixing conditions, but not in the strongly homogeneous Da^ <C 1 limit. 

Effects of particle inertia and gravitational settling on the droplet size distribution and vertical particle velocities 
have also been analysed in the DNS, with Stokes numbers not exceeding St v < 9 x 10~ 2 . The primary influence is on 
the initial evaluation of the negative tail of the Lagrangian supersaturation PDF and the resulting acceleration of the 
droplet evaporation. Within the parameter range studied, the mean droplet size was not modified by droplet inertia 
and settling. Likewise, mean droplet properties were not significantly altered with modest increases in Reynolds 
number. 

In this work, the whole study has been conducted in a small subvolume of clear air-cloud interface, specifically in 
a cubic box of dimensions up to L x = 51.2 cm. For such size, we expect the mixing to be dominantly homogeneous, 
with large-eddy Damkoehler numbers up to Da^— 2.7. This is exacerbated by the finding that the relevant phase 
relaxation time depends on the diluted droplet number density, so that the largest Dcil only reach 1.1. In the future, 
we intend to carry out a similar analysis in a larger box such that inhomogeneous mixing can take over at the larger 
scales of the flow. This will allow us also to make contact with recent large-eddy simulations^!. 

Furthermore, it can be expected that the inclusion of the active character of the temperature field will modify 
the droplet growth. In the present case the temperature was set constant thereby reducing to an advection-diffusion 
problem for the vapor field. Full thermodynamic consistency will require advection of temperature in the same 
turbulence, full consideration of latent heating associated with phase changes, and determination of the saturation 
vapor mixing ratio ratio as a function of the varying temperature field. Results of the present studies must therefore 
be interpreted in light of these simplifications, and considered to be a first step in building up to the full complexity 
of the cloud mixing problem. 
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